ASYMPTOTIC DYNAMICS OF 
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Abstract. We classify and predict the asymptotic dynamics of a class of swarming models. 
The model consists of a conservation equation in one dimension describing the movement of a pop- 
ulation density field. The velocity is found by convolving the density with a kernel describing 
attractive-repulsive social interactions. The kernel's first moment and its limiting behavior at the 
origin determine whether the population asymptotically spreads, contracts, or reaches steady-state. 
For the spreading case, the dynamics approach those of the porous medium equation. The widening, 
compactly-supported population has edges that behave like traveling waves whose speed, density and 
slope we calculate. For the contracting case, the dynamics of the cumulative density approach those 
of Burgers' equation. We derive an analytical upper bound for the finite blow-up time after which 
the solution forms one or more (5-functions. 
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1. Introduction. Biological aggregations such as fish schools, bird flocks, ungu- 
late herds, and insect swarms have drawn considerable attention from mathematical 
modelers in recent years. These animal groups - which for brevity we refer to simply 
as swarms - have implications for ecological dynamics, human food supply availabil- 
ity, disease transmission, and, on the longest spatiotemporal scales, evolution [T6|l20|. 
Increasingly, they serve as prototypes for the development of algorithms in robotics, 
engineering, and artificial intelligence [TH]. Furthermore, biological swarms are 
a rich and versatile source of pattern-forming behavior, taking on morphologies in- 
cluding vortices, advancing fronts, branched dendritic structures, and more exotic 
patterns [3[I7j. 

The emergent organization of swarms can be mediated by exogenous influences 
such as nutrients, light, or gravity, as well as by endogenous ones, namely social 
interactions between individuals. Since many species swarm even in the absence of 
meaningful external stimuli, one concludes that social interactions play a key role. The 
most important social forces are thought to be attraction, repulsion, and alignment [TJ 
|5J ^] . Attraction refers to the evolutionarily preprogrammed tendency of conspecific 
organisms to move towards each other, which offers benefits such as protection and 
mate choice, while repulsion refers to the tendency to move away, for instance, for 
collision avoidance |16j . Attraction and repulsion are driven by the relative locations 
of organisms. In contrast, alignment refers to the tendency of an organism to match 
the speed and orientation of its neighbors. 

As highlighted in the particular combination(s) of attraction, repulsion, and 
alignment that are included in a model strongly affect the classes of solutions observed. 
For example, [SJ [TD] have elucidated the importance of alignment in giving rise to 
diverse and exotic swarming patterns including pulses, breathers, and ripples. In 
contrast, models including attraction and repulsion as the only social forces have 
a history of several decades and an extensive literature, much of which is reviewed 
in [15]. These models typically give rise to groups that spread, contract, or reach 
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equilibrium |15 [ I23| , [22]. If organisms are self-driven in addition, milling and migrating 
groups may form [5J [T3J |T3| . 

In mathematical descriptions of swarms, a common modeling assumption is that 
social interactions take place in a pairwise manner, and that the effect of multiple 
organisms on a given organism can be determined via a superposition. Consider 
a swarm with a sufficiently large population such that the group is well-described 
by a continuum density p{x,t), as in [5] [TH [531 122 ^nd many others. Under the 
aforementioned modeling assumptions, social forces involve a convolution term of the 
form 



Here fs is a kernel describing the social influence of the population at location y on 
that at location x. Not only does the choice of social forces included in a model 
play a key role (as mentioned above), but the particular shape of the social kernel 
fs used to model a given social force can have a crucial affect on the dynamics of 
the group. For instance, the particular shape of the attractive-repulsive kernel used 
in [5J |S] determined whether groups collapsed into a dense group, formed a well- 
spaced vortex-like swarm, formed a ring-like structure, or took one of several other 
morphologies. 

If modelers are without explicit biological measurements giving an idea of a par- 
ticular organism's social kernel, they face a crucial question: in order to construct a 
model that gives the qualitatively correct swarming behavior, what kernel should be 
chosen? One might think that since the kernel is a function, it determines an infinite- 
dimensional parameter space, and so selecting a particular point in that space for 
one's model might be challenging. In practice, modelers typically choose a functional 
form that is presumed to be phenomenologically appropriate, for instance, a kernel 
fs that is exponentially decaying in space and has the correct near-field and far-field 
behavior. For a few examples, see Table 1 in [TS]. Even with such constraints, models 
may contain many parameters. For instance, there are five parameters controlling the 
social interactions used in [13], and at least eleven in [10] . 

In this paper we analyze a given class of swarming models with a general social 
interaction kernel, and we classify and predict the possible asymptotic dynamics. The 
class of models we consider is 



This equation describes a conserved continuum density field p{x,t) on the real line. 
The velocity v{x,t) depends exclusively on social interactions by means of a convo- 
lution with a kernel fs describing attraction and repulsion. In this paper, we focus 
solely on attractive-repulsive interactions, and hence do not consider social forces 
with an intermediate neutral zone as in, e.g., |11| : nor do we consider alignment. This 
model is kinematic, as opposed to dynamic, in which case the velocity would obey 
a momentum equation. As reviewed in |10j . social forces take place when animals 
communicate, either directly by auditory, visual, olfactory, or tactile senses, or indi- 
rectly, as mediated by chemical, vibrational, or other sorts of signals. A given type of 
communication may be unidirectional, as with visual sensing, or omnidirectional, as 
with auditory and olfactory sensing. Many organisms can process a combination of 




(1.1) 



Pt + {pv)x = 



(1.2a) 




(1.2b) 
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different input signals, which results in communication that is effectively omnidirec- 
tional |101ll8j . For this reason, in our one-dimensional model we choose fs to describe 
antisymmetric social forces, that is, we assume that fg is an odd function to ensure 
that distinct organisms exert equal and opposite forces on each other. Within the 



framework of (1.2), when sgn(x)/i(x) < then the effective social force is attractive, 



and when sgn(a;)/s(a;) > it is repulsive. Swarming models of the form (1.2) have 



been studied in [2J [3113] for specific choices of fs, and in a two-dimensional setting in 

A common choice for fs used in e.g. [HI UHl HI] and quite a few other studies is 
the Morse interaction force 

fs{x) = sgn(x) T-Fe-I"!/^ + e-l"lj . (1.3) 

Here, the first exponentially decaying term represents attraction with strength F > 
and characteristic length scale L > 0. The second term, of opposite sign, describes 
repulsion. The problem has been nondimensionalized so that the repulsive strength 
and length scale are unity. Figure [T3|a) shows a schematic example of (1.3) for the 



case F < 1, L > 1. The Morse function is, in fact, a member of the more general 
class of kernels 

fs{x) = sgn(x) [-Fg{\x\/L) + g{\x\)] , (1.4) 

where we scale the length and magnitude of g such that it has first moment equal 
to 2 and g{0~^) — 1. Here, g{x) is some suitable function: it could be a Gaussian, a 
compactly supported function, or one of many other choices. We analyze both the 
Morse function ( [1.3| a nd the more general class (1.4) in this paper. However, our 



goal is to analyze ( 1.2 1 with as few assumptions on fg as possible, so we also consider 



cases more general than (1.3 1 and (1.4) 



We have already assumed fs is odd. We make three additional, relatively weak 
assumptions in order to facilitate our analysis. First, fs has a finite first moment. 
This assumption is consistent with the idea that organisms should not interact at very 
long length scales because their range of sensing is limited. Second, fs is continuous 
and piecewise differentiable everywhere except for a finite jump discontinuity at the 
origin. The biological intuition that supports this assumption is as follows: for a given 
organism, the effect of other organisms in the far-field should vary continuously with 
distance. Small changes in distance should induce small, continuous changes in the 
social force. However, since fs is odd, it is discontinuous at a; = unless fs{x) — > 
'AS X I 0. This is a degenerate case which we exclude here since we expect organisms 
in close proximity to have nonzero effects on each other. Note that this assumption 
implicitly excludes the case of so-called "hard-core" forces that blow up at a; = 
|15j . Third, fs crosses for at most one value of We concentrate on the most 
biologically relevant case, when organisms are repelled at short distances (avoiding 
collision) and attracted at longer ranges (creating a tendency to form a swarm) . This 
means that when two organisms are within sensing range of each other, they have a 
unique pairwise equilibrium distance. For completeness, we will also consider other 
cases captured within our modeling framework, namely some cases where organisms 
only repel (i.e. fs>0 for x > 0) or only attract (i.e. < for x > 0), and briefly the 
unbiological case where there is attraction at short distances and repulsion at large 
distances. 



Our main results are as follows. Eq. (1.2), with fs as described above, has three 



possible asymptotic behaviors. The population density profile can spread, blow up, 
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Fig. 1.1. Schematic depiction of fs, the social interaction kernel in i f j . . (a) The Morse-type 
social force given by for the case F < 1, L > 1. (b) The piecewise-linear social force given by 

\2.2l\ for the case _F > 0. 



or reach steady state. Via long-wave and short-wave analyses, we predict how the 
asymptotic dynamics depend on fs- Specifically, the long-time behavior depends on 
two parameters which may be directly computed from fs'- the first moment and the 



limiting behavior at the origin. We perform numerical simulations of (1.2) to con- 
firm these predictions for two example kernels. For the spreading case, the dynamics 
approach those of the porous medium equation. The widening, compactly-supported 
population has edges that behave like traveling waves whose speed, density and slope 
we calculate. For the contracting case, the dynamics of the cumulative density ap- 
proach those of Burgers' equation. We derive an analytical upper bound for the finite 
blow-up time after which the solution forms one or more (5-functions. The case of 
steady-state solutions is studied in [T]. 

The remainder of this paper is organized as follows. In Section [2] we perform the 



long- and short-wave analyses of (1.2 1 to derive conditions on fs for the three pos- 
sible asymptotic behaviors, and we confirm these predictions with numerical simula- 
tions. We also derive (local) equations describing the asymptotic dynamics. Section 
|3] presents the spreading case in further detail, including an analysis of the traveling- 
wave-like behavior of the edge of the spreading group. Section |4] studies the blow-up 
case in more detail, including the analytical calculation bounding the finite blow-up 
time of the solution. We conclude in Section [Sj At the end of this paper, there are 
two appendices. Appendix |A] demonstrates conservation of mass and center of mass 



for (1.2 1. Appendix B gives an overview of a particle-based numerical method we 



developed to simulate the model. 

2. Asymptotic behavior of solutions. To demonstrate possible asymptotic 



behaviors of (1.2 1, we conduct numerical simulations using the Morse-type social in- 
teraction (1.3 1 as an example. Our simulation takes place on an infinite domain and 
uses a particle-based numerical method we have developed, described in Appendix [B| 
Simulations reveal three asymptotic behaviors, namely spreading, steady-state, and 



blow-up, as depicted in Figure 2.1 Figure |2.1[ a) shows a spreading solution, corre- 
sponding to a population that disperses to infinity. The population density profiles 
are compactly supported, with a jump discontinuity at the edge. The profiles ap- 



pear to be self-similar; we discuss this issue further in Sections 2.2 and [3] Figure 
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|2.1[ b) shows a steady-state solution, corresponding to a localized aggregation of the 
population. Again, the population density drops discontinuously at the edge of the 
support. Steady states of (1.2 1 are analyzed in [I]. Figures [2.1[c,d) show two cases 



of solutions where the density blows up, corresponding to populations with finite at- 
traction at short distances. In the first case, the solution forms a single clump. In the 
second case, the solutions form multiple, mutually-repelling clumps. These clumps 
are, in fact, (^-functions, as we discuss in Section |4] Our goal for the remainder of the 
present section is to derive conditions on a general social force fs to produce each of 
the aforementioned behaviors. To do this, we examine separately the long-wave and 
short-wave behavior of the system. 




(b) 



0.4 
0.3 
0.2 
0.1 





Fig. 2.1. Population density profiles governed by Ijl.S^ with Morse-type social interactions 
. . Arrows indicate the evolution of a profile over time. The asymptotic behavior of the model 
depends crucially on the choice of F, the relative strength of social attraction to social repulsion, 
and L, their relative characteristic length scales, (a) F = 0.2, L = 2. The compactly supported 
population eventually spreads to infinity, (b) F = 0.4, L = A. The population reaches a compactly- 
supported steady state, (c) F = 2, L = 2. The density profile blows up into a single clump, (d) 
F = 2, L = 0.5. The density profile blows up by forming mutually repulsive clumps. 



2.1. Long- wave behavior. We first consider the evolution of initial conditions 
that are wide and slowly-varying. Specifically, assume that p is initially long-wave, 
meaning p is localized near wave number fc = 0. We show that such initial conditions 
evolve, at least for a short time, according to the porous medium equation. To begin, 
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let us define the Fourier transform of a function h{x) as 

h{k) = T{h} = / h{x)e-''''' dx. 



(2.1) 



We apply the Convolution Theorem to (1.2b I to write the Fourier transform of the 
velocity as 



v{k)^3^{p*fs}^mfs{k). 

We next write fs{k) as a Taylor series 



(2.2) 



(2.3) 



Then, we express the n}^ derivative of /s at fc = in terms of the moments of fs ■ The 
n^^ moment of fs is 



z"/.(z) dz. 



Then, 



A(")(0) 



(-*)"Af„[/, 



Substituting (2.5el into (2.3 1, we obtain 



Z ^ 77 I 



C50 



n=0 



Since fs is antisymmetric, the even moments of fs vanish and we have 



n— ^ ^ 



or in physical space 



^ (2n+l)! (9a;2"+i 



(2.4) 



dfc" 7 


/,(x)e-^'=^ 

— OO 


fc=0 


(2.5a) 




fc=0 


(2.5b) 


(-i)" 


/ x"/,(x)e^*'=^dx 


(2.5c) 




^ — OO 




fe=0 


(-.)•■ j 


/•oo 

' a;"/,(x)dx 

— OO 




(2.5d) 



(2.5e) 

(2.6a) 
(2.6b) 

(2.7) 

(2.8) 
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If p varies on a lengthscalc much longer than the characteristic lengthscale of fs, its 
successive derivatives will be smaller and smaller. Assuming a nonzero first moment, 
we find that 



V{x) « -Mi[/s]/9:r + 0{pxxx)- 



(2.9) 



The successively smaller error terms in (2.8 1 correspond to higher-order (positive 
or negative) diffusion. With the velocity in this form, the governing equation ( |1.2[ ) 
becomes 



Pt = K(p2), 



(2.10) 



For K > 0, this is the well-known porous-medium equation. For certain initial 
conditions, a class of similarity solutions known as Barenblatt solutions are given by 



p{x,t) 



1 - 



X — Xq 



(2.11) 



-I + 



where we use the notation [u]+ = max{0,M}, and where M is the mass (c/. Ap- 
pendix A) , Xq is its center of mass and to is a parameter depending on the initial 
condition [25l|26]. Additionally, all initial conditions for (2.10) will approach this 
particular class of solutions asymptotically as t —^ oo. The solutions spread and grow 
wider without bound. For k < 0, (2.10) describes backwards diffusion; mathemati- 
cally, this problem is ill-posed. 



The case k > is asymptotically consistent; that is, long- wave states in (1.2 1 
will spread and therefore remain long-wave when k > 0. Eq. (2.10) will become an 



increasingly valid approximation of (1.2), justifying a posteriori the longwave expan- 



sion. However, when k < 0, long-wave initial conditions will contract until they can 
no longer be considered long- wave, at which point the approximations used above 
become invalid. Finally, we note that if k = (that is, if the first moment of fs 
vanishes), the above analysis does not hold, and we must retain higher-order terms 
in ( |2.7| in order to predict asymptotic behavior. 

For comparison purposes, we note that the root-mean-square (RMS) width of the 
Barenblatt solution (2.11) can be computed by first computing the second moment 



around the center of mass. 



Q 



which yields 



{x - Xo)^p{x,t) dy 



32/3 



34/3 



M5/3«;2/3(i + i^)2/3 



RMS = ^/qJm - " M^l^K^I^(t to)'/'. 
v5 



(2.12) 



(2.13) 



We verify this prediction in Section |3] below. 

2.2. Short-wave behavior. We now consider the evolution of initial conditions 
that are narrow and sharply-varying. We show that the cumulative density behaves, 
at least for a short time, according to Burgers' equation [4 . We first define the 
cumulative mass function: 



'(/'(a;, t) — p{z, t) dz. 



(2.14) 
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Note that since p> 0, ip{x) increases monotonically from a value of to M. We now 
use ( 1.2 1 to write 



ipt{x,t) = / Pt{z,t) dz 



{p{z,t)v{z,t))zdz 



= —p{x, t)v{x, t) 
= -il:x{x,t)v{x,t). 

That is, the cumulative mass function obeys 

i>t + vil:^: = 0. 



(2.15a) 

(2.15b) 

(2.15c) 
(2.15d) 



(2.16) 



To proceed, recall from Section [T] our assumption about the social interaction 
force fs, namely that fa is continuous and piecewise differentiable everywhere except 
for a jump discontinuity of size 2(3 at the origin. Following [4J, we write 

fs{x)^2pH{x) + g{x) (2.17) 

where (3 ^ Q, H is the Heaviside function, and g{x) is continuous and differentiable. 



Substituting (2.171 into (2.161 and using the fact that H * p = tp yields 



V'i + {2P^P + 5 * p) V:. = 0. 



(2.18) 



For convenience, and without loss of generality, let the (conserved) center of mass of 
p be at a; = 0. Assume that p is initially short-wave, so that p M in Fourier space. 
In this case, the term {g * p)il)x ~ Mgtp^. Since ^/ja; « p is short-wave, we may further 
approximate this term as Mg{0)tjja:- Using the fact that g{0) = (3 from (2.171 and 



substituting into (2.18), we have (approximately) that 



V-t + (2/3V + M/3) V. = 



(2.19) 



for short-wave solutions. This is Burger's equation with an additional constant veloc- 
ity term. This term may be eliminated by a simple change of variables, for instance 
letting — AI/2 to obtain 



^At + 2(3^^p, = 0, 



P = lim/s(a:). 



(2.20) 



We now invoke standard results for Burgers' equation |24-. Since ip is monotonically 
increasing in x, ip will contract and form a shock when /3 < and will spread when 
(3 > 0. Moreover, because V-'x = Pi ^ shock in ^ is manifest as a (5-function in p, 
and so we expect that p will blow up when (3 < and spread when /3 > 0. In fact, 
under mild conditions on fs , 4 shows global existence of solutions for (3 > and give 
examples of finite-time shock formation for (3 < 0. 

The case /? < is asymptotically consistent; that is, short-wave initial conditions 
in (1.2 1 will contract and therefore remain short wave when /3 < 0. Eq. (2.201 will 



become an increasingly valid approximation of (1.2 1. However, when /? > 0, short- 



wave initial conditions will spread until they can no longer be considered narrow, at 
which point the approximations used above will fail to hold. Finally, we note that 
if /3 = 0, we are in the degenerate case where fs is continuous at the origin. In this 



case, the leading order approximation for (1.2 1 would involve antiderivatives of the 
cumulative mass function -ip. 
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2.3. Predicting qualitative behavior. From the results in Sections |2.1| and pT2| 

we expect short waves to blow up when /3 < and spread when /? > 0. Similarly, we 
expect long waves to contract when k < and spread when k > 0. When short waves 
blow up, we expect the short-wave instability to override the long-wave behavior. 
Thus, there are three possible cases, one of which has two sub-cases. We summarize 
these below. 

(A) When [3 > and k > 0, both long and short waves expand, leading to spread- 
ing solutions of the type shown in Figure 2.1 'a). The asymptotic dynamics of 
the density are governed by the porous medium equation (2.101. We analyze 
this case further in Section [3l 

(B) When (3 > and k < 0, short waves spread while long waves contract, leading 
to steady-state solutions of the type shown in Figure [2J^ b). We analyze these 
solutions in depth in yy. 

(C) When /3 < 0, short waves contract and solutions blow-up regardless of the 
value of K, leading to solutions of the types shown in Figure 2.1 'cd). The 
asymptotic dynamics of the cumulative density are governed by Burgers' 
equation (2.201. We analyze this case further in Section |4] 

Note that, plausibly, there could be other asymptotic behaviors that we have not 
discovered. However, the value of f3 at the origin governs whether organisms are 
repulsive or attractive at short distances and clearly governs the formation of clumps. 
If /3 > 0, p must spread to at least a width where the long-range attractive forces 
play a significant role. Our interpretation of k is as a measure of whether long-range 
attraction can balance short range repulsion (the case when k < 0) or if the short- 
range repulsion always dominates (when k > 0). While more exotic behaviors might 
be possible with more exotic choices of fs - say with multiple bands of attraction and 
repulsion - the classification above captures the behaviors observed with the simple, 
biological choices of fs considered in this paper. 

As an example we consider the class of social forces ( 1.4 1, for which k = 1 — FL^ 
and (3 = 1 — F. (Note that the regime F > 1, L > 1 corresponds to purely attractive 
social forces, and the regime F < 1, L < 1 corresponds to purely repulsive social 
forces.) We expect to see blow-up when F > 1, spreading when F < 1/L^ and F < 1, 
and steady-state solutions when 1 > F > These predictions are indicated in 

Figure pT2| which shows F-L parameter space. The blow-up boundary (3 — 1 ~ F — 
is the solid line and the spreading/ steady-state boundary k = 1 — FL^ = is the solid 
curve. The symbols in Figure |2.2| summarize the results of numerical simulations 



conducted for the particular case when fs is the Morse function ( 1.3 1. The theoretical 



curves divide the numerical results as expected. Our Figure [2T2] is similar to "phase 
diagrams" showing the linear stability and statistical mechanical H-stability of other 
swarming models with Morse-type social forces in [5J |5] . 

In the blow-up regime, we in fact observe two different types of blow-up in the 
numerical simulations; the boundary between these is indicated as the broken vertical 
line. The particular form of the blow-up depends on the long-range character of the 
social force fs- For the class of kernels (1.4 1, when L > 1, fs{x) ^ 0~ as a; ^ oo 
so social forces are attractive at long distances. In this case, the entire mass of 
the system eventually collapses into a single (5-function. In the other case i < 1, 
fs{x) ^ 0+ as X ^ oo and social forces are repulsive at long distances (a behavior 
which does not have an immediate biological interpretation). Blow-up still occurs due 
to the contraction of short waves. However, the long-range repulsion means that the 
solution does not aggregate into one clump. Instead, multiple (5-functions form which 
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repel each other and move apart. 
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Fig. 2.2. Different dynamical regimes of the governing equation l[1.2y in F-L parameter space 
for social forces of the form F is the relative strength of attraction to repulsion, and L is the 

relative length scale. The solid horizontal line at F = 1 indicates (3 = in ^2.20}^ and marks the 
theoretical boundary above which solutions blow up. The solid curve F = indicates k = in 

\2. 1 Q[ ) and marks the theoretical boundary between spreading and steady state solutions for F < 1. 
Results of numerical simulations using Morse-type social forces are indicated by the symbols: 

spreading (o), steady-state (x), and blow-up (Q ando). The (partial) line L = 1 indicates the critical 
case separating whether the long-range behavior is attractive or repulsive. In the blow-up regime with 
long-range attraction (L> 1), the solution forms a single S-function (o). With long-range repulsion 
(L < I), multiple, mutually-repelling 5-functions form (O). 



To verify our analytical results further, we consider a second example with a social 



force not of the form of ( 1 .4 1 , namely 



fs{x) = sgn(x) 



1^ 

l-L 

.0 



'x\ 



F < |a;| < L 
Z < |a;| < 1 
\x\ > 1 



1^ 

l-L 



(2.21) 



where F G M. and L S (0, 1). (Note that the regime F < corresponds to purely 
attractive social forces within the range of sensing.) A schematic picture is shown in 
Figure [Ll|b) for the case F > 0. For a; > 0, this function consists of the compactly 
supported, piecewise linear function passing through the points (0,F), (L,— 1), and 
(1,0). For X < the function is the odd extension. The parameter F plays a role 
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somewhat similar to in ( 1.3 1 in that it determines whether the kernel is attractive 

The parameter L plays a role somewhat similar to 
Eq. (lOT) differs from Ol 



or repulsive for short distances. 

L in that it determines a characteristic length scale, 
in that the kernel is compactly supported rather than decaying, is linear rather than 
exponential, and by construction is attractive (negative) at intermediate distances 
regardles s of p arameter choices. 

For (2.211, K — {FL"^ — L — l)/6 and (3 ~ F. We expect to see blow-up when 

F 




< JD, spreading when F>L +L >0, and steadv-st ate s olutions when 



L . These predictions are indicated in Figure 



2.3 



which is simi- 



< F < L' 

lar to Figure |2.2| except that now we use a social force given by (|2.2ip rather than 



( 1.3 ). As before, numerical simulations produce spreading solutions and steady-states, 
both with jump discontinuities at the edges, as well as solutions that blow up. The 
theoretical predictions for these different regimes (curves) again divide the numerical 
results, as expected. For this example, only single ^-function blow-up occurs because 
/s < at intermediate distances and our initial conditions have sufficiently narrow 
support. Since fg is compactly supported, initial conditions that are sufficiently wide 
(or consist of sufficiently distant, separated groups) would blow-up into multiple 6- 
functions that are stationary, rather than mutually repelling. 



3. Spreading Solutions. When k > in (2.101 and /3 > in (2.201, solutions 



will spread. As discussed in Section [2] the density profile grows wider, the long-wave 
approximation (2.10) will become increasingly accurate, and so we expect solutions 
to approach Barenblatt's solution ( |2.11| . Figures 3.1 and 3.2 confirm this prediction. 



Figure 3.1 compares Barenblatt's solution (2.11 1 to numerical simulations of ( 1.2 1 



0.2 > in (2.101 



using the Morse-type social force (|l.3|) with F — 0.2 and L = 2. For these parameters 

K 



The broken line represents Barenblatt's solution (2.111 for a 
density profile with unit mass. Under the similarity transformation 



p{x) = —p{jx), 7 = maxp(a;) 

7 X 



(3.1) 



the spreading Barenblatt profile collapses to a single curve. The solid lines represent 



snapshots from the numerical simulation of (1.2 1. We apply to these numerical snap- 



shots the same rescaling (3.1 1. We take the initial condition to be a rectangular pulse 



with unit mass. The direction of increasing time is indicated by the arrow in the fig- 
ure. As time increases, the numerical profiles, as expected, approach the Barenblatt 
profile. We explore this approach further in Figure (3.2), which compares the root- 



mean-square (RMS) width of the solution. From (2.13), the RMS width should grow 
as t^/^. As predicted, the RMS width of the numerical solution (circles) approach the 
theoretical curve (line) on the log-log plot. 

Apart from the shape of the solution, we also wish to study the jump disconti- 
nuities at the edge of the spreading swarm. At the left edge of the swarm, we might 
expect a spreading solution to behave locally like a fixed wave profile traveling to 
the left (and similarly at the right edge). Therefore, we seek a traveling- wave solu- 
tion to (1.2). At the left endpoint, we look for a traveling- wave solution of the form 
p{x,t) — g{x + ct), where g{z) = for all z < 0. Under these assumptions, (1.2 1 



reduces to 



(3.2) 
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Fig. 2.3. Differ ent dy namical regimes of the governing equation 1.2) in F- L para meter space 
for the social force 2.21). The horizontal line at F = indicates (3 = in \ 2.20) and marks 
the theoretical boundary below which solutions blow up. The solid curve F = + indicates 
K = in ^2. 1 Q[ ) and marks the theoretical boundary between spreading and steady state solutions. 
Results of numerical simulations are indicated by the symbols: spreading (o), steady-state (x), and 
single-S -function blow-up (o). 



Integrating both sides of this equation, 

{c + v)g = 0. 



(3.3) 



The constant of integration is zero because the left-hand side vanishes for negative z. 
Hence, wherever g is nonzero, ~c — v. That is. 



g{i)fs{z-z)dz forz>0. (3.4) 
We proceed with a quantitative analysis for the example case of Morse-type social 



interactions (1.3 1. Writing out explicitly and taking derivatives with respect to z 
on both sides (which eventually facilitates transformation of the integral equation into 
an ODE) yields, after some rearranging. 



{F~l)g{z) 



5(5) 



F 



-\z-z\IL 



-\z-z\ 



dz 



(3.5) 



To ensure that the exponential terms are linearly independent, we assume F ^ Q 
and L ^ 1. Then, to solve this integral equation, we apply the differential operators 
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Fig. 3.1. Spreading solutions to using the Morse-type social interaction . 3[ ) with F = 0.2 

and L = 2. We choose as an initial condition a rectangular pulse with unit mass. Snapshots of 
the evolving profile are rescaled according to the similarity transformation ^3. These evolving 
solutions are the solid curves and the arrow indicates the time evol ution . As predicted, the numer- 
ical solutions approach the idealized Barenblatt similarity solution {g. J which has been similarly 
rescaled and is shown as the broken curve. 



^1 — dzz — 1 and £2 ~ L'^dzz — 1 to both sides. The left-hand side becomes 

= (F - l)L^g""{z) + (1 - F + L2 - FL^)g"{z) + {F ~ l)g{z), 
and the right-hand side becomes 



g{z)CiC2 



10 

/"OO 



F 



^-\z-z\IL _ g-|z-z| 



dz 



{L^-F)g"{z) + {F-\)g{z). 



Hence, the integral equation (3.4) reduces to the ODE 

g""{z)-a^g'\z)^Q 

where 



a 



L^{l-F)' 



(3.6a) 
(3.6b) 



(3.7a) 



g{z) ■ (-2) • [{F - L^)S"{z - z) -I- (1 - F)5{z - z)] dl (3.7b) 



(3.7c) 
(3.8) 

(3.9) 




Fig. 3.2. RMS width of the solutions in Figure 
solution and has equation RMS = 3^/''5~^/^(t + io)^ 
represent the numerical solutions to 
rate. 



The solid line represen ts Ba renblatt's 
0.544tl/^ at large times \2.13 ). Circles 
1.^ which asymptotically approach the Barenblatt spreading 



Because we are in the spreading regime by assumption, 1 — FL? = k > and 1 — F = 
/? > 0, so > and thus the coefhcient on g"{z) is strictly negative. Integrating the 
ODE twice yields 



g"{z)~a^[g{z)-Az-B] = Q 



which has general solution 



-l{z) ^Az + B + Ce-"^ + De° 



(3.10) 



(3.11) 



The traveling wave cannot grow exponentially as z — > oo, as this would imply expo- 
nentially growing mass flux (which is proportional to the product of the speed and 
the derivative of the profile) as the wave translates to the left, so we choose D — 0. 
To find A, B, and C we plug (3.11 ) into (3.4 1 and simplify to obtain 



A- 



B ■ 



2{FL^ - 1) - FL^e-^/^ 

FLe-'/^ 
FL 



(3.12) 



e 



1 - aL 



1 

( 

1 — a 
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Since {1, e" 



^} are linearly independent, we can match like terms and solve 



the resulting three algebraic equations for A, B, and C to obtain 
A = cAo, Aq 
B = cBo, Bo 
C — cCq, C'o - 



1 


2(1 






1 


2(1 


-FL2) 




1 


2(1 





L + l 



{aL- 1) 1 



(3.13a) 
(3.13b) 
(3.13c) 



which determines a traveling-wave solution for each wave speed c. 



Figure 3.3 shows an example of the traveling left edge of the spreading swarm 
studied in Figures |3.1| and |3.2[ We plot three snapshots of the numerically spreading 



solution and superpose the analytical solution given by (3.111 and (3.13). The two 
are in good agreement close to the edge of the support where the traveling wave 
calculation above is expected to be valid. 




Fig. 3.3. Evolving left edge of the spreading soluti on stu died in Fig ures \3.1\ and \3.^ We su- 
perpose the analytical solution (dashed curves) given by \3.11}^ and _? 3[ ) on the numerical solution 
(solid curves). The two are in good agreement close to the edge of the support where the traveling 
wave calculation in Section\^is expected to be valid. 



To verify (3.11) and (3.13) further, we predict the relationship between the in- 
stantaneous speed of a traveling front, the size of the jump at the edge, and the slope 
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of the density at the edge. In particular, the size of the jump is 

g(0) = c(Bo + Co), (3.14) 

and the slope at the edge is 

g'{0)^A-aCe~"'\^^^^c{Ao~c^Co). (3.15) 

A similar analysis holds at the opposite edge of the swarm. We expect these relations 
to hold only for large t since the solution must be sufficiently wide and slowly varying 
near the endpoints for it to locally approximate a traveling wave. For several values 
of F and L, we tested these predictions by tracking the speed, jump in density, and 
slope at the endpoints over time. Figure |3.4| shows an example that confirms the 
traveling wave predictions. This example corresponds to the same spreading profile 



studied in Figures |3.1| and |3.2[ Denote the location of the left edge of the swarm by 

(3.16) 



Xe- We plot three ratios involving quantities computed at the edge, namely 



c(Bo + Co)' c(Ao-aCo)' p,(xe)(Bo + Co) 
where we take as the values of p{xe) and Px{xe) their limit approaching from the inside 



of the support. Each of the three quantities in (3.161 approaches unity as t ^ oo as 



predicted by (3.141 and (3.151 



4. Contracting solutions. We now consider the case when (3 < and solu- 
tions blow up due to short-wave contraction. Biologically, this means that at short 
distances, organisms are attracted to each other leading to clumping. As discussed 
in Section [2] the density profile p blows up by forming one or more 5-functions, or 
equivalently by the cumulative density ip forming shocks. The space-time plot Figure 



4.1 



shows an example. Lines represent contours of '0 in the simulation of (1.2 1. The 
value is coded by shading, indicating the characteristics of the hyperbolic problem. 
As we expect, the characteristics intersect and form a shock after suflicient time, cor- 



responding to blow-up of p. The authors of [5] study (1.2) for the case when /s < 
and rigorously show blow-up in finite time. The blow-up profiles for certain other 
instances of fs are studied in [4], which also shows finite-time blow-up for /3 < 0. In 
Section |2] we showed blow-up for /? < 0, regardless of initial conditions, for any fs 
satisfying our prior assumptions 

When the initial condition is a single, sufficiently narrow pulse, we can approxi- 
mately predict when the solution will form a 5-function. Let a(t) and b{t) denote the 
position at time t of the left and right edges, respectively. Then, for a{t) < z < b{t) 
and b{t) — a{t) sufficiently small, note that 

fs{a{t)-z)> min /,(r) (4.1) 

a(t)-b{t)<r<Q 

= max Mr) (4.2) 

a<r<b{t)-a{t) 

min \fs(r)\. (4.3) 

Here, we know fs{b — a) < because b — a is small and positive, and because 
f3 = linir,. iQ fs{x) < 0, with continuous except at the origin (c/. Section [2]). For 
convenience, define 



q{t)= min |/.(r)|. (4.4) 

0<r<b{t)-a{t) 
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10"' 10° io' lo" 10^ 10* 



t 



Fig. 3.4. Edge behavior of a spreading solution. This example corresponds to the profile studied 
in Figures \S.l\ and \3.t^ The squares, triangles, and circles indicate, respectively, the three ratios in 
j^g. j 6[ ) . These quantities all approach unity (horizontal dotted line) for large t, which indicates that 
the asymptotic dynam ics o f the e ndpoints obey the results of the traveling wave analysis. Specifically, 
as predicted by Ii3.14^ and \3.15^ , the jump in density and the slope of the density profile at the edge 
are both proportional (via known constants) to the speed at which the edge moves. 



We find a bound for the velocity at the left endpoint: 



a'{t) = v{a{t),t) 



p{z,t) fs{a{t) - z) dz 



a(t) 



Mt) 

>q{t) / p{z,t)dz 

Ja(t) 

= Mq{t) 
> 0. 

A similar argument shows that the velocity at the right endpoint satisfies 

b'{t) < -Mq{t) < 0. 



(4.5) 
(4.6) 

(4.7) 

(4.8) 
(4.9) 

(4.10) 



Thus, the endpoints approach each other. From (4.4 1 it follows that q{t) must be 
non-decreasing, and consequently the endpoints must be accelerating towards each 
other, or at least moving towards each other at a constant velocity. 
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0.2 0.4 0.6 0.8 1 1.2 1^1 



t 



Fig. 4.1. Space-time plot of the cumulative mass ■(/> showing blow-up of the density p under the 
dynamics of We use the Morse function J j . 3^ with attractive strength F = 2 and attractive 

length scale 2, for which f3 = —1 in ^2.2C^ . Contours of ^ appear as lines, and the value is coded 
by shading, indicated the characteristics of the hyperbolic problem. As we expect, the characteristics 
intersect and form a shock, which means that p blows-up by forming a S-function as predicted. 



Let t* denote the time at which all the mass of the system enters a single S- 
function. We can find an upper bound for t* by noting 

b'{t) < 6'(0) < -Mq{0), a'{t) > a'(0) > Mq{0). (4.11) 

Since the endpoints are initially separated by a distance &(0) — a(0) and each is moving 
towards the other at a minimum speed Mq{0), this gives the bound 

bjO) ajO) 

- 2Mq{0) • ^^-^^^ 

Furthermore, just before the solution forms a (5-function, by similar argumentation, 
the velocities of the endpoints will be M\f3\ at the left endpoint, and — M|/3| at the 
right endpoint. If /3 = but attraction dominates at small distances, a careful analysis 
of q{t) suggests that either finite-time or infinite-time blow-up can occur. 



We have studied blow-up in numerical simulations of (1.2 1; results appear in 



Figure 4.2 We use the Morse-type social force (1.3 1 with L — 2. Blow-up time is 
shown as a function of F for two different initial conditions, namely a rectangular 
pulse of width 0.1 (circles) and one of width 0.2 (squares). Both sets of data closely 
match the analytical upper bound, indicated as solid and broken curves respectively. 
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Fig. 4.2. Blow-up times t* for with Morse-type interactions ^77^. We set L = 2 and 

vary F . Data correspond to two sets of initial conditions, namely a rectangular pulse of width 0.1 
(circles) and one of width 0.2 (squares). The actual blow-up times are well-approximated by the 
analytical upper bound 0773), indicated as the solid and broken curves. 



5. Conclusions. In this paper we have studied the swarming-type equation ( 1.2 1 
with the goal of predicting how the asymptotic dynamics depend on the social inter- 
action force /s. We consider the class of social interactions that are antisymmetric 
(describing isotropic interactions), have finite first moments, and are continuous ex- 
cept at the origin. From a long-wave and a short-wave analysis, we showed that two 
parameters computed directly from fa determine the asymptotic dynamics, namely 
the first moment 2k, and the limit approaching the origin from the right (3. When 
/? > and K > 0, long and short waves both expand and as i — > oo, the dynamics 
approach those of the porous medium equation. The shape of the profile and its 
spreading rate approach those of Barenblatt's well-known solution. For the case of 
Morse-type interactions, we calculated a quantitative relationship between the edge 
speed, edge density, and edge slope. An interesting question for analysts is to deter- 
mine conditions on fg such that this self-similar solution is a global attractor. It is 
clear that /? > and k > are necessary, and we conjecture that our restriction that 
fs has at most a single zero crossing (for a; > 0) is sufficient. 

In the case /3 > 0, k < 0, long waves contract and short waves expand. In this case, 
the system asymptotically reaches a steady state whose shape is highly nontrivial. 
In [I], we analyze these solutions in detail both for the original governing equation ( 1.2 1 
as well as for the case when (1.2b I contains an additional term describing exogenous 
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forces acting on the population (for instance, the effects of gravity, hght, or a nutrient 
field). 

When (3 < 0, short waves contract. Regardless of the long-wave behavior, solu- 
tions blow up as i ^ cx). In this limit, the dynamics of the cumulative density obey 
Burger's equation and form shocks; hence, the density forms one or more (5-functions. 

There are several clear directions for future work. First, with the definition of 
the social force fg correspondingly modified, many results from Section [2] could be 
extended to the case of higher dimensions. In particular, the predictions of the bound- 
aries of the different dynamical regimes would be of interest. Second, though we have 
studied in [T] the effect of exogenous forces on steady-state solutions, we have not 
studied their effect on the spreading and contracting solutions that are the focus of 
the present paper, nor their potential effect on shifting the dynamical regime bound- 
aries. Other extensions would include the addition of alignment forces and loosening 
of the restriction of omnidirectional communication, both of which would require an 
appropriate reformulation of the governing equations of motion. Such a study might 
shed light on the role parameter choices play in the models investigated in [9l [10] . 

We hope that our present study will guide mathematicians, biologists and engi- 
neers who wish to construct swarming models with particular behaviors that either 
mimic those observed in nature or are desirable qualities for control of robotic or vir- 
tual agents. Our results suggest that although many functional forms can be imagined 
for the social forces, only a few types of qualitative behavior manifest for this class 
of model. From another perspective, selecting a particular functional form to model 
social interactions is less important than choosing the parameters in that model to 
manifest a desired behavior. Finally, this study suggests that there are a number 
of characteristics inherent to kinematic models, namely spreading parabolic profiles, 
small-scale clumps, and groups with jump discontinuities at the edge, which may be 
used to diagnose when this class of models is appropriate. 
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perceptive comments which greatly improved this paper. 

Appendix A. Basic properties of the continuum equation. 



We demonstrate that the governing equation (1.2 1 conserves mass and center of 



mass. We take the social force fg in (1.2b) to be odd, which describes the case of 



isotropic social interactions. That is, organisms in disparate locations have a social 
effect on each other that is equal in magnitude and opposite in direction. 



Conservation of mass follows from the fact that (1.2 1 is formulated as a conser 



vation law. To see this explicitly, define the mass of the system at time t as 

/oo 
p{x,t)dx, (A.l) 
-oo 

and so 

dM 



dt 



Pt dx = - / (pv)^ dx^ - [pv] = 0, (A.2) 

-oo J —oo 



assuming the density decays to zero as x — > ±cx). We denote the mass of the system 
at any time as M = M(0). 
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Intuition suggests that the center of the mass of the system should also remain 
fixed because the antisymmetry of the social interaction force fg- To verify this, we 
consider the center of mass at time t: 

1 f°° 

= y xp{x,t)dx. (A.3) 

Then, 



dx 1 

'dt ~ M 

1 



xpt dx (A. 4) 



x{pv)xdx (A. 5) 



L f J x=-oo 

Assuming the density vanishes at ±oo 
dx 1 



dt M 
1 

~ M 
1 

^ M 



pv dx (A. 7) 

/oo 

p{y,t)fs{x-y)dydx (A.8) 
-OO 

p{x,t)p{y,t)fs{x~y)dydx. (A.9) 



-oo 

oo /"OO 



oo — oo 



Relabeling the variables of integration and invoking the antisymmetry of the social 
interaction force. 



dx _ I 
'dt ~ M 
_ 1 



— OO J — OO 
oo pOO 



p{y,t)p{x,t)fs{y - x) dxdy (A. 10) 

p{x,t)p{y,t)fs{x ~ y)dydx. (A.U) 



oo ./ — oo 



Hence, dx/dt — —dx/dt, and thus dx/dt — 0. That is, the center of mass is stationary. 
Appendix B. Numerical method. 

Our numerical solution of ( |1.2| hinges on a correspondence with a discrete model 
that approximates it. Ref. j3j shows that a discrete model of the type we will derive 
converges to the continuous model under fairly weak assumptions. 

Our correspondence works as follows. Consider a continuous distribution pc{x,t) 
with total mass M. For ease of notation, we suppress time dependence for the re- 
mainder of this paragraph. Define the cumulative density function 

M^)= f Pc{s)ds (B.l) 

where the dummy coordinate xq is taken to the left of the support of pc- We seek a 
discrete approximation of N (5-function point-masses each of mass m = M/N. That 
is, 

N 

Pd{x) ^^m6{x - Xi). (B.2) 
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The associated cumulative density function ij^d is 

'O X < xi 

^.(x) = <'™[^/^ + ^'"^)^ (B.3) 
im Xi < X < Xi+i, i — 1, . . . , N — I 

X > xn 

where we have used the convention that integrating up to a (5-function yields half the 
mass of integrating through it. To establish a correspondence between the discrete 
and continuum problems, we require that ipcixi) = ^d{xi), which in turn determines 
the point-mass positions Xi. As N oo for fixed M, this step function ipd converges 
uniformly to ipc- The correspondence goes in the opposite direction as well. If we 
begin with an ensemble of (5-functions pd, we can find the corresponding cumulative 
density tpd, interpolate to approximate ipc, and differentiate to find an approximate 
continuous density pc- 

With this correspondence established, we now describe our numerical method. 
Given an initial condition pc{x,0), we determine the corresponding discrete density 
Pd{x,0) and the initial point-mass positions Xi(0) as described above. Substituting 



(B.2| into the governing equations (1.2) yields a system of N ordinary differential 



equations 

^ = ^m,/,(x,(t)-x,(i)). (B.4) 
See |15| for an introduction to discrete swarming models of this type. We then solve the 



differential equations (B.4 1 numerically. From the new point-mass positions Xi{t), we 
then reconstruct 'ipd{x), ipc{x), and pdx, t), again using the correspondence described 
in the preceding paragraph. 

Our numerical scheme has three sources of error. First, there is error associated 



with the integration of the ordinary differential equations (B.4 1. This error is easily 
controlled. We use the Matlab routine ODE45 for the numerical solution. The second 
source of error is interpolation error in the construction of pc from the location of the 
point masses. We interpolate "0^ to construct tpc and then differentiate the polynomial 
analytically to obtain pc- The error is 0{N~^). The third source of error comes 



from the approximation of the integral in the velocity term (1.2b I. We perform this 
quadrature using the point-masses for collocation, with an error of 0{N~^). For a 
full description of the numerical method and further details of the error analysis, see 

m. 
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